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Abstract. - More than 35 million high-energy muons collected with the MACRO detector at 
the underground Gran Sasso Laboratory have been used to search for flux variations of different 
nature. Two kinds of studies were carried out: a search for the occurrence of clusters of events and 
a search for periodic variations. Different analysis methods, including the Scan Statistics test and 
the Lomb-Scargle spectral analysis have been applied to the data. 



Introduction. — The high energy muon events collected by the MACRO apparatus at 
the average depth of 3800 m.w.e. represent one of the most extensive records of such kind of 
data. These data can be used to search for time variations of periodic and stochastic nature, 
as it was done extensively by using arrival times of Extensive Air Shower [1]. Variations in 
the underground muon flux may be due to different causes of galactic, solar and terrestrial 
origin. The common problem for this type of searches is to determine whether an observed 
effect has occurred by chance or if it signals a departure from a pure random distribution. 

Bursts of underground muon events may be originated by violent events as GRB on the 
time scale of few seconds or by sudden metereological variations on the scale of several hours. 
Periodic modulations are related to the cooling of the upper atmosphere during the night 
(winter) producing daily (seasonal) variations of the atmospheric density and therefore in 
the underground muon flux. A sidereal modulation may be introduced by the Earth roto- 
traslation in a isotropical cosmic ray distribution (Compton-Getting effect). 

MACRO was a multipurpose modular apparatus with 6 supermodules equipped with 
liquid scintillators, limited streamer tubes and nuclear track detectors [2]. The complete 
detector was a nearly "closed box" with a total length of 76.7 m; it had a lower part 4.8 
m high, with 10 horizontal layers of streamer tubes, two layers of scintillators and seven 
layers of rock absorber. The upper part ("attico") was 4.5 m high and contained 4 layers of 
streamer tubes, one layer of scintillators and the electronics. The streamer tube system was 
used for tracking, and provided two independent views: the wire view and the strip view. 
The latter employed 3 cm wide aluminium strips at 26.5° with respect to the wire view. 

MACRO studied atmospheric neutrinos and their oscillations [3,4], various aspects of 
cosmic ray physics and astrophysics [5], searched for GUT Magnetic Monopoles [6] and 
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Fig. 1: Distribution of the muon rate per run, in units of /i/hour. We accepted runs within 3rr from 
the average of the gaussian fit superimposed to the data. The small peak on the left corresponds 
to runs with the detector partially in acquisition (1/3 of the supermodules were off). The tail 
at the left of the primary peak (~1% of the total) corresponds to runs where one module was in 
maintenance and the data taking was still progressing for particle searches. 

other exotica [7]. Some interruptions of different kinds occurred during data taking, either 
randomly (e.g. power outages), or regularly (e.g. maintenance), so appropriate statisti- 
cal methods have to be applied and particular care should be used in choosing periods of 
stationary conditions. 

In the following we discuss the results of the searches for periodic variations and for time 
clustering of muon events. 

Burst search. — Using the MACRO data we performed a new analysis on muon time 
series in order to find isolated clusters of events. The MACRO Collaboration published 
the results of a previous analysis obtained with another method, namely the study of the 
time interval distribution [8]. For each muon arriving at time to the distribution of the 
time interval elapsed between the first muon to and the next five muons: t^-to, i=l,...,5 was 
calculated. 

For the present analysis we considered data recorded by the streamer tube system in the 
time interval November 1991-May 2000 and we selected the data with the following criteria: 

- run duration longer than 1 hour; 

- streamer tube efficiency larger than 90% for the wire view and 85% for the strip view. 
The streamer tube efficiency was obtained using the sub-sample of muon tracks crossing all 
the 10 lower streamer tube planes; 

- all 6 super-modules in acquisition. 

- acquisition dead time smaller than 2.5% for the whole detector. 

- runs having muon counting rates that deviated more than 3cr from the average were 
removed (Fig. 1). It was checked that this cut did not bias the analysis: for the Scan 
Statistics, the expected burst durations are small compared to the run durations; for the 
periodicity search, the expected modulation amplitudes are negligible compared to the fluc- 
tuations around the average number of muons per hour. 

- we required a single reconstructed track in each one of the two projective views. This 
cut excludes muon bundles from the analysis and provides a clean sample of single muon 
events. 

The total number of runs surviving the cuts was 6113 corresponding to 3.5 • 10 7 muon 
events. This sample was used to search for bursts of cosmic ray particles coming from the 
whole upper hemisphere. In this analysis, we looked for bursts of events using the so-called 
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"scan statistics" method. This is a bin-free method and it provides unbiased results (see [9] 
and references therein). 

Let us consider an interval A = [A, B] of a continuous variable x and a Poisson process 
with density A. We call scan statistics (SS) the largest number of events found in any 
subintcrval of [A, B] of length lu 

S{w) ee max {Y x (w)} (1) 

A<x<B— uj 

where Y x (w) denotes the number of events in the subinterval [x,x + w]. The probabil- 
ity P that a statistical fluctuation would produce a burst of events as large as k can be 
approximated by [9]: 

P = 1 - Q*(k; ij)L, 1/L) ~ 1 - Q* [Q* 3 /Q* 2 ] L - 2 (2) 

where tp = Au>, L = A/w and 

Q* 2 = [F(k-l,i,)f -(k~l)p{k 1 4>)p{k~2,^) 

-Qe-l-il>)p(k,il>)F(k-3 t il>) (3) 
Ql = [F(k-1,^)] 3 -A 1 +A 2 + A s -A 4 (4) 



with 



Ai = 2p(k,xp)F(k-l,ip){(k-l)F(k-2,ip) 

-i()F{k-3,ip)} 
A 2 = 0.5[p(fc,^)] 2 {(fc- l)(fc-2)F(fc-3,V) 

A 3 = ^p(2fc-r,V)[^(r-l,^)] 2 

r=l 
fe-1 



A 4 = ^2p(2k-r,ip)p(r,ip){(r-l)F(r~2,'i(j) 
-^F(r-3»} 

In the above formulas F(k, if)) denotes the cumulative distribution 

F(A,^)=5^p(i J V) ; p(i,i>) = e-^ (5) 

i=0 

and F(k, tp) = for k < 0. 

We used the SS in the following way: for each run i, let L4i,£?i] be the time interval 
ranging from the start to the end of the run. The interval was scanned counting the number 
of muons inside a "time window" w of fixed length. Let ki be the maximum number of events 
recorded during the scan. For each run, according to Eq. 2, we computed the probability 
Pi that a statistical fluctuation would produce a burst of events as large as ki. The only 
parameter to be fixed "a priori" is the scanning window w, whose size must be chosen on 
the basis of astrophysical considerations (e.g. time duration of typical HE burst events). We 
tried different sizes (w = 30 s, 5 min and 15 min) and for each the probability distribution 
Pi (i=l, N rura ) was analyzed. 

Fig. 2 shows the S(w) distributions for the three different values of the scanning window 

w. 

In Fig. 3 we show the probability distribution for the 6113 runs surviving the analysis 
cuts: w — 30 s (Fig. 3 top), 5 min (Fig. 3 middle) and 15 min (Fig. 3 bottom). No 
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Fig. 2: Distribution of the S(w) SS variable for the three sliding windows used in this work. Each 
entry of the histograms corresponds in the horizonthal scale to the maximum number of events 
falling within a given time window for each run. 

significant deviations from the null hypothesis was found. The analyses of unusual runs 
with probabilities smaller than 5 ■ 10~ 4 have shown that "bumps of events" were located 
near the beginning or the end of the runs. The structures at Logio(Prob) > —0.4 in the 
plot with w=30 s are related to the time window's width comparable with the muon rate 
of ~0.2 muons/s. In this case non-Poissonian sources due to background, dead time, etc. 
start to dominate and cause a deviation from the expected distribution. 

Periodicity search and spectrum analyses. — The Fourier amplitude spectrum 
analysis is a powerful technique that allows a blind search for regular/persistent fluctuations 
in time series [10]. Such a technique, however, requires the input data to be sampled at 
evenly spaced intervals; data gaps of variable length and occurring randomly in the series 
produce spurious contributions to the power that can mimic the presence of a periodicity. 

The Lomb-Scargle method [11] mitigates this effect, even in the case of very long data 
series. Moreover, as indicated in Ref. [12] , it allows to evaluate the significance of the "peaks" 
(signal) with respect to the null hypothesis. The data used in this analysis were collected 
with the full detector and with strict selection criteria. The muon events were binned in 
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Fig. 3: SS probability distributions for the selected runs. In the upper panel a time window w — 30 
s was used; w — 5 min in the central panel and w = 15 min in the bottom panel. 
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Fig. 4: Lomb power p as a function of the frequency [days -1 ] for experimental data (upper panel). 
Note the high peak at ~2.7-10~ 3 corresponding to 365 days (seasonal flux variation). In the lower 
panel the results of a Monte Carlo simulation having the same noise level of the real data with 
seasonal, solar diurnal and sidereal waves [13, 14] added. 



15 min time intervals and bins deviating by more than 3er from the monthly average rate 
were discarded. The total number of time-bins used was 160242 corresponding to 58% of 
the whole sample. 

The results of our analysis are shown in Fig. 4. We compare the spectrum obtained 
for the real data (upper panel) with a Monte Carlo simulation having the same noise level 
and time intervals distributed according to the sequence of the original series (lower panel). 
Seasonal, solar diurnal and sidereal waves were added in the simulated series with the 
amplitudes found in Refs. [13,14]. Real data have a power distribution which decreases 
with frequency -1 / 2 up to f=l day -1 : this behavior, mainly related to the dead time of the 
apparatus, was not included in the simulation. The most striking feature of the spectrum 
(note the logarithmic vertical scale) is the large peak at ~ 2.7 • 10~ 3 corresponding to the 
seasonal variation of the muon flux as found in Refs. [13, 14]. 

Fig. 5 shows a frequency region around the solar diurnal frequency where are indicated 




Fig. 5: The frequency region around the solar diurnal wave. The arrows mark its position and also 
the sidereal and anti-sidereal peaks. The upper panel refers to real data, the lower panel to Monte 
Carlo simulation. 
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the frequencies corresponding to the sidereal and anti-sidereal waves, both for data (upper 
panel) and Monte Carlo simulation (lower panel). In order to "clean" the power spectrum 
from unphysical frequencies, we performed a "sliding average" both for data and Monte 
Carlo: each time-bin was rescaled according to the formula: 

N' = Ni ~ ^ (6) 
N{At) { 1 

where iVj is the original bin content and N( At) is the average content in the time range ±Ar 
(we chose Ar = 1 day). The peak at frequency 1 day -1 has a statistical significance of ~2.3 
a; the statistical significance assuming an oscillatory behavior is ^3.4 a. We looked also in 
the frequency regions around the sidereal and anti-sidereal waves. A signal corresponding 
to the sidereal variation is observed, but peaks of similar size (or even larger) are present 
elsewhere in the spectrum. The claim that the sidereal and solar diurnal waves are real is 
based on their occurrence at a frequency of "a priori" interest and on the stability of its 
amplitude and phase with time: this is an indication of an independent observation of the 
same modulations found with the "folding" method described in Ref. [14]. The amplitudes 
and the probabilities for the null hypothesis computed with the folding method and with 
the Lomb-Scargle method are in fair agreement. 

Conclusions. — We analyzed the time series of MACRO muons using two complemen- 
tary approaches: search for bursts of muon events and search for periodicities in the muon 
time distribution. The SS method was used in the first case and the Lomb-Scargle method 
in the second case. The two techniques complete early analyses performed with "folding" 
methods in searching for periodicities and time differences for burst events. The seasonal 
modulation (Fig. 4) was confirmed. The signals at the positions of the solar diurnal and 
sidereal modulations are confirmed, even if with smaller significance (Fig. 5). No other 
deviations from the expected distributions were found. 
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